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SOUND  EXPOSURE  LEVEL  PREDICTION  FOR  IMPULSE 
SOUND  SOURCES  ABOVE  VARIABLE  TERRAIN 


1  INTRODUCTION 


Background 

Blast  noises  from  the  use  of  heavy  weapons  on  Army  training  ranges  may  propagate  over  several 
kilometers,  adversely  affecting  surrounding  communities.  Being  able  to  calculate  the  effects  of  Army-type 
blasting  would  help  Army  planners  to  locate  ranges  or  schedule  training  where  and  when  it  would  cause 
the  least  disturbance.  However,  the  influence  of  ground  surface  makes  it  difficult  to  calculate  the  sound 
wave  propagation  or  the  prediction  of  sound  exposure  level  (SEL)  of  a  point-acoustic  source.  Even  if  a 
rigid  ground  surface  is  assumed,  the  reflection  and  diffraction  of  sound  waves  by  realistic  terrain  still 
makes  the  calculation  very  complicated.  It  is  similarly  difficult  to  calculate  electromagnetic-wave 
propagation  near  the  ground  surface.  Luebbers  (1985)  overcame  this  difficulty  by  u^ing  a  propagation 
model  based  on  Geometrical  Theory,  or  ray  theory,  of  Diffraction  (GTD)  for  electr'magnetic  waves. 
Although  acoustic  waves  usually  have  much  lower  frequencies  than  elecromagnetic  waves,  the  ray  theory 
of  sound  waves  has  been  an  efficient  and  realistic  mcxiel  for  the  propagation  of  acoustic  waves.  Later 
research  extended  the  model  to  include  acoustic  wave  propagation  over  hills  and  ridges  (Houshmand 
1990).  Although  the  acoustic-wave  GTD  model  represents  a  significant  step  towards  more  accurate 
numerical  calculations,  several  important  features  must  be  added  to  the  model  before  it  can  be  applied  to 
realistic  situations. 

The  GTD  program  uses  a  two-dimensional  model.  The  terrain  is  modeled  as  a  piecewise,  linear, 
two-dimensional  structure  of  rectangular  plates.  (Although  the  terrain  is  only  two-dimensional,  the  source 
and  receiver  can  be  located  anywhere  above  the  surface.)  In  modeling  realistic  terrain,  the  surface  can 
be  specified  as  an  arbitrary  assembly  of  plates,  each  having  any  aspect  and  slope. 


Objectives 

The  objectives  of  this  work  were  to  extend  the  GTD  propagation  model  to  include  realistic  terrain, 
to  use  the  improved  model  to  man  'ue  spatial  SEL  prediction  of  a  given  explosive  source  above  a 
specified  terrain,  and  compare  those  predictions  to  measured  data. 


Approach 

A  literature  search  was  conducted  to  locate  an  appropriate  mathematical  strategy  for  modeling 
sound-wave  propagation  over  realistic  terrain.  The  (GTD)  model  of  Luebbers  was  adapted  to  include 
realistic  terrain,  broadband  sounds  (blast  waves),  and  frequency  weighting  factors.  The  improved  model 
was  interfaced  with  the  Geographic  Resources  Analysis  Support  System  (GRASS)  to  calculate  the  SEL 
from  impulse  sources,  given  the  map  elevation  contours. 


Mode  of  Technology  Transfer 

The  results  of  this  study  will  be  used  as  baseline  guidance  in  future  work  in  predicting  propagation 
over  rough  or  variable  terrain. 
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2  THE  RAY  THEORY 


Since  the  discovery  of  Newton’s  three  laws  of  motion  in  1687,  people  have  been  accustomed  to 
coasider  motion  as  something  that  occurs  along  a  certain  path.  The  propagation  of  waves  in  an  extended 
medium  has  been  considered  as  motion  of  perturbation  along  a  path  of  rays;  this  is  known  as  ray  theory. 
Under  the  ray  theory,  concepts  such  as  refraction,  reflection,  and  diffraction  are  easily  understood.  The 
path  of  rays  can  be  mechanically  traced  by  Snell's  law. 

In  geometrical  theory,  the  rays  along  which  the  waves  propagate  in  a  homogeneous  medium  are 
rectilinear  and.  upon  striking  the  boundary  between  two  such  media,  are  reflected  back  into  the  original 
medium  or  refracted  into  another  medium  according  to  certain  simple  rules.  The  rules  are  valid  in  the 
approximation  only  if  path  differences  of  the  order  of  a  wavelength  can  be  neglected. 

From  the  standpoint  of  wave  theory,  the  concept  of  a  wavefront  plays  a  central  role  in  the 
geometrical  theory  or  ray  theory.  A  wavefront  is  any  moving  surface  along  which  a  waveform  feature 
is  being  simultaneously  received.  It  is  not  necessarily  assumed  that  the  amplitude  along  a  wavefront  is 
coastant,  or  that  the  wavefront  is  planar.  However,  geometrical  theory  assumes  that  the  amplitude  varies 
only  slightly  over  distances  comparable  to  a  wavelength  and  that  the  radii  of  curvature  of  the  wavefront 
are  substantially  larger  than  a  wavelength. 

Although  audible  sound  waves  have  much  larger  wavelengths  (17  cm  to  17  m)  man,  for  example, 
the  wavelength  of  visible  light,  which  occupies  a  range  in  the  neighborhood  of  50  pm,  there  is  still  wide 
use  of  geometrical  acoustics  or  ray  acoustics,  which  developed  as  an  important  branch  of  acoustical  theory 
(Pierce,  1989).  Ray  acoustics  is  an  approximation  in  the  same  sense  as  is  geometrical  optics,  in  that, 
within  the  range  of  a  wavelength,  both  the  medium  and  the  boundary  are  homogeneous.  The  influences 
of  the  actual  spatial  fluctuations  of  the  medium  or  of  the  boundary  are  counted  as  noise  or  perturbation 
on  the  results. 

Luebbers  (1985)  developed  a  computer  model,  the  “GTD  propagation  model,”  to  compute  the 
electric  field  from  a  dipole  antenna,  the  source  point,  radiating  above  a  piecewise  linear  two-dimensional 
terrain  structure  where  each  linear  segment  of  the  structure  may  have  arbitrary  permittivity,  conductivity 
and  surface  roughness.  The  model  computes  the  total  complex  electric  field  at  the  field  point  (the 
observation  point)  by  locating  and  summing  all  of  the  rays  of  the  electromagnetic  waves  that  can 
propagate  from  the  source  dipole  to  the  field  point  either  directly  or  by  interaction  with  the  linearized 
terrain  profile.  The  individual  rays  contributing  to  the  electric  field  may  be  obtained  separately,  so  that 
the  user  can  obtain  information  as  to  how  most  of  the  energy  is  reaching  the  field  point,  for  example,  by 
reflection  from  a  certain  area  or  diffraction  over  a  particular  ridge. 

The  Luebbers’s  ray  propagation  model  was  intended  for  use  in  the  frequency  range  of  100  MHz  to 
10  GHz.  Houshmand,  ct  al.  (1990)  extended  the  GTD  propagation  model  to  acoustic  wave  propagation. 
A  vector  electric  field  calculated  is  replaced  by  a  scalar  sound  pressure.  However,  Houshmand’s  model 
accounts  only  for  a  restricted  model  for  the  terrain,  such  that  the  terrain  is  composed  of  an  assembly  of 
parallel  ridges  or  valleys.  This  structure  can  be  specified  as  height  as  a  function  of  one  coordinate  only: 
that  running  perpendicular  to  the  ridges.  This  report  relaxes  the  assumption  to  a  more  realistic  model. 
In  this  new  model,  the  terrain  can  be  specified  as  an  arbitrary  function  of  position  limited  only  by  the  grid 
size.  The  SEL  computation  may  still  be  accomplished  with  the  theory  of  ray  acoustics.  Any  inhomo¬ 
geneity  of  the  terrain  within  the  order  of  the  grid  size  is  considered  as  “roughness"  of  the  terrain  surface. 
Although  ray  acoustics  may  be  not  as  good  in  this  case  as  wave  theory  because  of  the  low  frequency  of 
sound  waves,  the  geometrical  method  should  not  be  ignored  as  a  way  to  estimate  SEL. 
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Two-Dimensional  GTD  Propagation  Model 

Figure  1  shows  an  example  of  a  parallel  terrain  profile  that  has  already  been  linearized  for 
conducting  a  ray  calculation  (Luebbers,  1985).  In  this  example,  15  points  make  up  a  p.ofile  consisting 
of  14  connecting  segments.  Remember  that  in  any  direction  x,  perpendicular  to  the  paper,  the  terrain  is 
infinite.  Each  point  represents  an  edge,  and  each  segment  represents  a  plate  of  terrain  surface.  After 
inputting  the  terrain  parameters,  the  location  of  the  wave  source,  and  observation  points,  and  the 
frequency,  the  model  traces  all  the  rays  that  can  exist  and  judges  those  rays  that  can  propagate  from  the 
source  point  to  the  observation  point.  The  complex  amplitudes  of  ail  the  rays  that  can  be  received  at  the 
observation  point  are  summed  up  to  obtain  the  final  wave  amplitude  and  phase.  Several  of  the  types  of 
rays  included  in  the  model  are  illustrated  in  Figures  2  to  5.  These  figures  show  that  the  direct  ray  of 
Figure  2  propagates  directly  from  the  source  to  the  observation  point,  the  reflected  ray  of  Figure  3  is 
reflected  once  from  the  terrain,  and  that  the  reflected-diffracted  ray  of  Figure  4  is  reflected  once  and 
subsequently  diffracted  once  on  its  path  to  the  observation  point.  In  Figure  5,  we  see  a  reflected- 
diffracted-refiected  ray  that  reflected  from  Plate  K  to  Edge  e,  then  diffracted  from  Edge  c  to  Plate  M,  and 
finally  reflected  from  Plate  M  to  the  observation  point,  a  human  listener,  for  example.  The  GTD  model 
considers  16  types  of  rays  (Table  1). 

The  theoretical  basis  for  determining  direction  of  reflection  of  rays  is  simply  to  make  the  angle  of 
incidence  equal  to  the  angle  of  reflection.  The  amplitude  and  phase  of  the  reflected  electromagnetic  waves 
is  determined  by  the  complex-valued  reflection  coefficient,  computed  using  the  angle  of  incidence, 
electrical  constants  of  the  reflecting  segment,  and  a  roughness  factor  representing  height  variability  in  local 


Figure  1.  Sample  Terrain  Profile  That  Allows  All  16  Types  of  Rays  Included  in  Table  1  of  the 
GTD  Model. 
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Figure  4.  The  Reflected-Diffracted  Ray  Reflected  From  Plate  K  to  Edge  e,  and  Diffracted  From 
Edge  e  to  the  Observation  Point. 


Figure  5.  The  Reflected- Diffracted -Reflected  Ray  Reflected  From  Plate  K  to  Edge  e.  Diffracted 
From  Edge  e  to  Plate  M,  and  Reflected  From  Plate  M  to  the  Observation  Point. 
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Table  1 


Types  of  Rays  Included  in  the  GTD  Model 


JDIR - Direct  ray; 

JREF - Reflected  rays,- 

JRR - Reflect  twice; 

JRD - Reflect,  then  diffract; 

JRRD - Reflect  twice,  then  diffract; 

JRDR - Reflect,  diffract,  then  reflect; 

JDIF - Diffract  once; 

JDR - Diffract,  then  reflect; 

JDRD - Diffract,  reflect,  then  diffract; 

JDD - Diffract  twice; 

JDDR - Diffract  twice,  then  reflect; 

JDRR - Diffract,  then  reflect  twice; 

JRDD - Reflect,  then  diffract  twice; 

JRRR - Reflect  three  times; 

JDDD - Diffract  three  times; 

JRDDD - Reflect,  then  diffract  three  times. 


terrain,  such  as  variability  caused  by  vegetation,  uneven  ground,  or  waves  (if  propagation  is  over  water). 
The  roughness  factor  is  used  to  modify  the  reflection  coefficient  to  account  for  imperfect  reflection  caused 
by  local  terrain  roughness.  Specifically,  the  modified  reflection  coefficient  for  a  rough  surface  used  by 
the  GTD  model  is  defined  as  the  plane-wave  reflection  coefficient  for  specular  reflection  from  a  flat 
smooth  surface  multiplied  by  a  surface  roughness  factor  5.  The  surface  roughness  factor  can  be  found 
from  the  standard  deviation  AG  of  the  terrain  elevation  along  each  piecewise-Iinear  section  of  terrain 
defining  the  profile  as  follows: 


with 


AO  = 


:  e  -M2 

[Eq 

4kAG  .  ,, 

x  s,n7 

[Eq 

U 

2] 


where  y  is  the  grazing  angle  and  A,  is  the  wavelength.  In  Eqs  1  and  2,  AO  represents  the  phase  shift 
between  the  shortest  and  longest  reflected  path  from  an  ensemble  of  rays  that  strike  the  rough  surface. 

For  diffraction,  the  basic  approach  taken  is  to  approximate  each  change  in  slope  in  the  piecewise- 
Iinear  approximation  to  the  terrain  profile  as  a  wedge.  Wedge  diffraction  solutions  have  been  available 
since  the  beginning  of  the  century,  but  have  not  been  widely  used  for  propagation  path  loss  predictions 
due  to  their  complexity.  A  numerically  simple  and  elegant  solution  to  wedge  diffraction  which  serves  as 
a  canonical  problem  in  the  geometrical  theory  of  diffraction  was  published  by  Keller  (1962)  30  years  ago. 
However,  it  was  not  suitable  for  general  application  to  radio  propagation  over  diffracting  edges,  since  it 
fails  in  the  vicinity  of  shadow  boundaries,  i.c.,  when  the  source,  edge  and  field  point  lie  near  a  straight 
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line.  This  shortcoming  was  eliminated  by  Kouyoumjian  and  Pathak  (1974)  in  their  formulation  of  a 
wedge  diffraction  coefficient,  which  is  valid  at  shadow  and  reflection  boundaries. 


As  originally  formulated,  the  diffraction  coefficients  were  limited  to  perfectly  conducting  wedges. 
A  heuristic  extension  that  allows  the  approximate  treatment  of  finitely  conducting  locally-rough  wedge 
diffraction  has  been  made  by  Rojas-Teran  and  Burnside  (1981),  Chamberlin  and  Luebbers  (1982)  and 
Luebbers,  et  al.  (1982).  Luebbers  (1985)  further  made  modifications  to  the  usual  wedge  diffraction 
coefficients  nessessary  to  include  finite  conductivity  and  local  surface  roughness.  The  amount  of 
diffracted  energy  re-radiated  from  an  edge  is  determined  by  the  diffraction  coefficient,  which  is  a  function 
of  the  wedge  angle  and  the  incident  and  diffracted  ray  geometries.  Note  that  all  these  considerations  and 
modifications  fit  both  electromagnetic  waves  and  acoustic  waves. 


Ray  Selection  Algorithms  for  Two-Dimensional  Situations 

A  computer  is  needed  not  only  to  be  able  to  evaluate  individual  ray  amplitudes,  but  also  to  check 
for  blockage  and  to  search  for  all  possible  combinations  of  each  type  of  ray  interactions  with  a  given 
terrain  profile  and  then  to  combine  them  to  get  the  total  field  at  the  observation  point. 

Direct  Ray 

The  direct  ray  is  the  most  significant,  if  it  exists,  and  is  defined  as  the  ray  that  travels  directly  from 
the  original  source  to  the  observation  point.  This  ray  can  exist  only  if  there  is  no  terrain  high  enough  to 
block  its  path  of  travel.  Since  the  actual  terrain  is  approximated  by  many  finite  flat  plates,  the  computer 
program  must  logically  test  the  ray  for  blockage  by  any  of  these  plates.  If  the  ray  is  not  blocked,  then 
its  associated  wave  field  is  calculated,  and  it  becomes  the  first  wave  field  contribution  included  in  the  total 
wave  observed. 

Singly-Reflected  Ray 

A  singly-reflected  (R)  ray  is  defined  as  a  ray  reflected  once  from  a  plate,  where  the  incident  angle 
of  the  ray  with  respect  to  the  plate  is  equal  to  the  reflected  angle.  To  test  for  the  existence  of  the  ray,  the 
first  step  is  to  find  the  image  source  location.  This  image  source  is  uniquely  determined  by  the  original 
source  and  the  unit  normal  to  the  plate.  After  the  image  source  is  found,  the  ray  is  assumed  to  radiate 
from  this  source  to  the  observation  point,  with  its  amplitude  multiplied  by  the  appropriate  reflection 
coefficient.  The  second  step  is  to  determine  whether  or  not  the  ray  has  an  interception  point  within  the 
physical  boundaries  of  the  plate.  Also,  the  ray  must  not  be  blocked  by  the  other  plates.  Once  this  ray 
passes  all  of  the  tests,  its  wave  field  is  calculated  and  vectorially  added  to  previous  contributions.  This 
computer  algorithm  is  again  performed  on  the  next  plate  until  all  of  the  plates  have  been  included.  The 
final  result  is  the  total  singly-reflected  wave  field. 

Singly-Diffracted  Ray 

Diffraction  is  the  phenomenon  that  occurs  when  the  wave  ray  is  incident  on  any  discontinuity  in  the 
geometry  of  an  object,  i.e.,  a  sharp  comer  or  an  edge.  One  of  the  important  stages  in  the  process  of 
finding  the  diffracted  wave  field  is  to  locate  the  diffraction  point,  Qe  (Figure  6).  Qe  is  the  point  where 
the  angle  between  the  line  from  the  source  to  Qe  and  the  diffracting  edge  is  equal  to  the  angle  between 
the  diffracting  edge  and  the  line  from  Qg  to  the  observation  point.  In  other  words.  Qg  is  the  point  that 
satisfies  the  Law  of  Edge  Diffraction  (Fermat’s  Principle). 
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Figure  6.  The  Singly-Diffracted  Ray  Diffracted  From  Edge  e  to  the  Observation  Point. 


Once  Qe  is  found,  the  logical  test  for  ray  existence  is  performed.  Let  Plate  L  and  Plate  L+l  form 
a  wedge  e,  as  shown  in  Figure  6.  Edge  e  is  tested  for  the  existence  of  a  diffraction  point.  The  computer 
program  separates  the  terrain  into  two  regions.  The  First  region  includes  the  First  plate  to  plate  L-l;  the 
second  region  is  from  plate  L+2  to  the  last  plate.  In  the  first  region,  the  computer  program  tests  for 
blockage  using  the  source  location  and  Qg  as  the  two  end  points  of  a  ray.  In  the  second  region,  Qe  and 
the  obsevation  point  are  the  two  end  points,  and  similar  blockage  tests  are  performed.  If  the  ray  satisfies 
these  blockage  tests,  the  diffracted  ray  is  evaluated,  and  the  complex  amplitude  is  added  to  the  total.  If 
the  ray  was  blocked,  the  ray  is  abandoned.  The  program  repeats  the  procedure  until  all  of  the  edges  have 
been  tested  for  possible  diffracted  rays.  The  final  result  is  the  total  singly  diffracted  field. 

Doubly-Reflected  Ray 

The  doubly-reflected  (RR)  ray  is  one  that  is  reflected  from  two  different  plates.  An  example  of  this 
ray  is  shown  in  Figure  7,  where  the  ray  is  reflected  from  Plate  L  to  Plate  M,  and  reflected  off  Plate  M 
to  the  observation  point.  The  existence  determination  algorithm  starts  with  the  location  of  the  doubly- 
imaged  source.  First,  the  image  of  the  original  source  is  found  with  respect  to  Plate  L,  treat  this  image 
source  as  though  it  were  the  actual  source,  and  find  its  reflection  with  respect  to  Plate  M.  The  RR  ray 
radiates  from  this  doubly-imaged  source  to  the  observation  point,  but  with  its  amplitude  multiplied  by  the 
two  reflection  coefficients.  A  line  from  the  doubly-imaged  source  through  the  observations  point  intersects 
Plate  M  at  the  reflection  point  (CM).  If  CM  is  located  outside  the  physical  boundaries  of  Plate  M,  the 
RR  ray  does  not  exist.  But  if  CM  lies  within  the  boundaries,  the  reflection  point  on  Plate  L  (CL)  must 
be  found.  The  first  imaged  source  and  CM  form  a  ray,  whose  reflection  point  CL  must  lie  on  the  Plate 
L.  If  CL  and  CM  exist,  the  blockage  tests  are  performed.  In  summary,  for  the  RR  ray  to  exist,  inter¬ 
ception  points  CL  and  CM  must  lie  on  both  plates,  and  the  lines  of  sight  from  the  original  source  to  CL, 
from  CL  to  CM,  and  from  CM  to  the  observation  point  must  not  be  blocked.  If  these  conditions  are  met, 
the  RR  ray  wave  amplitude  and  phase  are  calculated  Mid  added  to  the  previous  contributions.  This  process 
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is  repeated  until  all  of  the  possible  plate  combinations  have  been  tested.  The  final  result  is  the  total 
doubly-reflected  field. 

Reflected-Diffracted  Ray 

The  reflected-diffracted  (RD)  ray  is  the  ray  that  is  reflected  once  from  a  Plate  K  and  then  diffracted 
once  by  an  Edge  e.  An  example  of  such  a  ray  is  shown  in  Figure  8. 

The  existence  of  a  reflection  point  (CK)  and  the  diffraction  point  (Qe)  is  the  first  requirement  for 
an  RD  ray  to  exist.  There  are  at  least  two  possible  methods  to  determine  the  existence  of  these  two 
points.  The  first  method  involves  creating  images  of  the  edge  under  test  and  the  observation  point  with 
respect  to  an  extension  of  the  reflecting  plane  being  tested  (K  in  this  case).  The  geometries  involved  in 
the  first  method  are  shown  in  Figure  9a.  With  the  source  location,  image  edge,  and  image  observation 
point,  the  location  of  the  diffraction  point  on  the  (image)  edge  is  calculated.  Finally,  if  a  line  between 
the  source  and  the  image  of  Qc  intersects  the  Plane  K,  CK  exists,  and  the  ray  is  tested  for  blockage. 

The  second  method  (see  Figure  9b)  is  to  image  the  source  rather  than  the  diffraction  point.  For  this 
method,  a  line  between  the  source  image  and  Qg  intersecting  Plane  K  determines  the  existence  of  Qe  and 
CK.  Due  to  computational  considerations,  the  latter  method  is  considered  desirable  and  is  used  in  the 
model,  although  both  methods  give  identical  results. 

If  CK  and  Qe  exist,  the  blockage  algorithm  is  performed.  If  the  ray  is  not  blocked,  RD  ray  wave 
amplitude  and  phase  are  calculated  and  added  to  the  previous  contributions.  The  amplitude  is  also 
multiplied  by  the  appropriate  reflection  coefficient. 

This  process  is  repeated  until  all  plate  and  edge  combinations  for  the  RD  ray  have  been  tested.  The 
final  result  is  the  total  reflected-diffracted  field. 
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Figure  7.  The  Doubly-Reflected  Ray  Reflected  From  Plate  M,  and  Reflected  From  Plate  M  to 
the  Observation  Point. 
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Figure  8.  The  Reflected-Diffracted  Ray  Reflected  From  Plate  K  to  Edge  e,  and  Diffracted  From 
Edge  e  to  the  Observation  Point. 
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Two  Methods  for  Determining  the  Existence  of  a  Diffraction  Point  (Qe)  for  a  Possible 
Reflected-Diffracted  Ray  Reflected  From  Plate  K  and  Diffracted  at  Edge  e. 
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Diffracted-Reflected  Ray 

The  diffracted-reflected  (DR)  ray  is  diffracted  from  an  edge  and  reflected  (off  a  plate),  which  is  not 
one  of  the  two  plates  that  make  up  the  edge.  The  ray  geometry  is  illustrated  in  Figure  10.  The  first 
requirement  for  a  DR  ray  to  exist  is  that  the  diffraction  point  and  reflection  point  must  exist  There  are 
two  possible  methods  to  determine  the  existences  of  these  two  points. 

The  first  method  involves  reflecting  the  source  and  the  edge  (e  in  this  case)  with  respect  to  an 
extension  of  the  plate  under  test  (L  in  this  case).  Figure  1  la  depicts  the  geometry  involved  in  this  method. 
The  image  of  the  diffraction  point  is  determined  from  an  image  source  location,  an  image  Edge  e  and  the 
observation  point.  Using  the  diffraction  point  image  and  the  observation  point  as  two  end  points,  together 
with  a  plane  that  contains  Plate  L,  reflection  point  (CL)  is  calculated.  Finally,  Qg  is  determined  by 
reflecting  the  image  of  the  diffraction  point  with  respect  to  Plate  L. 

The  second  method  is  to  reflect  the  observation  point  with  reject  to  the  plane  containing  Plate  L 
(see  Figure  lib).  With  the  source  location,  the  edge,  and  the  image  of  the  observation  point,  Qe  is 
determined  using  the  same  logic  as  for  the  singly-diffracted  ray.  CL  is  found  by  using  Qe  and  the  image 
of  the  observation  point  as  two  end  points  of  a  line  that  intercepts  an  extension  of  Plate  L. 

A  comparison  of  these  two  methods  indicated  that  the  latter  method  required  fewer  computational 
steps,  hence  less  computation  time.  Therefore,  the  latter  method  is  used. 

After  Qe  and  CL  are  determined  to  exist,  the  second  requirement  is  that  the  line  of  sight  from  the 
source  location  to  Qe  from  Qe  to  CL,  and  from  CL  to  the  observation  point  must  not  be  blocked.  If  the 
line  of  sight  is  not  blocked,  the  DR  ray  wave  amplitude  and  phase  are  calculated  and  added  to  the 
previous  contributions.  The  DR  ray  is  also  multiplied  by  the  reflection  coefficient  for  Plate  L.  This 
process  is  repeated  until  all  edge  and  plate  combinations  for  the  DR  ray  have  been  tested.  The  final  result 
is  the  total  diffracted-reflected  field. 


Figure  10.  The  Diffracted- Reflected  Ray  Diffracted  From  Edge  e  to  Plate  L  and  Reflected  From 
Plate  L  to  the  Observation  Point. 
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Doubly-Diffracted  Ray 

The  doubly-diffracted  (DD)  ray  is  diffracted  by  two  different  edges.  An  example  of  this  ray  is 
shown  in  Figure  12,  where  the  first  edge  is  formed  by  connected  Plates  K  and  K-l,  and  the  second  edge 
is  formed  by  Plate  L  and  Plate  L-l. 

Before  starting  the  blockage  test,  the  two  diffraction  points  DK  and  DL  have  to  be  determined.  If 
both  DK  and  DL  exist,  testing  for  blockage  of  the  line  of  sight  from  the  source  to  DK,  DK  to  DL,  and 
DL  to  the  observation  point  is  performed.  If  there  is  not  line  of  sight  blockage,  the  DD  ray  exists,  and 
its  associated  wave  field  components  are  calculated.  This  process  is  repeated  until  all  of  the  edges  have 
been  tested.  The  final  result  is  the  total  doubly-diffracted  field. 

Triply-Reflected  Ray 

The  triply-reflected  (RRR)  ray  is  defined  as  any  ray  reflected  from  plates  three  times,  with  the 
incident  and  reflected  angles  on  each  plate  equal.  Figure  1 3  shows  an  example  of  such  a  ray. 

To  test  for  the  existence  of  the  ray,  the  first  step  is  to  find  the  image-source  location  sound  source 
with  respect  to  the  first  reflecting  plate  (K  in  this  case).  The  method  of  determining  the  image  source 
location  is  the  same  as  for  previous  rays  involving  reflection.  Once  the  image  location  is  found,  using 
it  as  a  source,  a  second  image  with  respect  to  Plate  L  is  located.  Then  using  the  second  image  as  if  it 
were  the  source,  a  third  image  point  with  respect  to  Plate  M  is  located. 

The  second  step  is  to  determine  whether  or  not  reflection  points  exist  within  the  physical  boundaries 
of  Plate  K,  L,  and  M.  This  can  be  ascertained  by  using  the  third  image  and  the  observation  points  as  two 
points  on  a  line.  This  line  must  intercept  Plate  M;  the  interception  point  (RM)  and  the  second  image 
location  are  then  used  as  two  points  on  a  line  that  must  intercept  Plate  L.  Again,  by  using  the  interception 
point  on  Plate  L  (RL)  and  the  first  image  location  as  two  end  points  on  a  line,  the  third  interception  point 
(RK)  is  found  and  must  be  located  on  Plate  K.  If  any  one  of  the  three  reflection  points,  RK,  RL,  and  RM 
is  not  located  within  the  appropriate  plate  boundary,  the  RRR  ray  does  not  exist  and  the  next  reflecting 
plate  M+l  is  considered. 

If  three  reflection  points  exist,  the  third  step  is  to  determine  whether  or  not  the  ray  that  travels  from 
the  source  to  RK,  from  RK  to  RL,  from  RL  to  RM,  and  from  RM  to  the  observation  point  is  blocked  by 
the  other  plates. 

If  the  RRR  ray  passes  all  blockage  tests,  its  wave  amplitude  and  phase  are  calculated  and  added  to 
the  previous  contributions.  This  algorithm  is  performed  on  all  combinations  of  plates.  The  final  result 
is  the  total  triply-reflected  field. 

Reflected-Reflected-Diffracted  Ray 

The  reflected-reflected-diffracted  (RRD)  ray  is  defined  as  a  ray  reflected  by  two  different  plates,  then 
diffracted  from  an  edge,  which  must  not  be  part  of  the  two  reflecting  plates,  to  the  observation  point  An 
example  of  this  ray  is  given  in  Figure  14. 

The  first  step  to  test  for  the  existence  of  the  ray  is  to  find  the  image  source  position  with  respect 
to  Plate  K.  Once  the  first  image  point  is  found,  the  second  image  point  must  be  calculated  using  the  first 
image  location  and  Plate  L.  Using  the  second  image  location  and  an  edge  (Edge  e  in  this  case)  the 
diffraction  point  is  determined. 
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Figure  11.  Two  Methods  for  Determining  the  Diffracted-Reflected  Ray  Path. 
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Figure  12.  The  Doubly-Diffracted  Ray  Diffracted  From  Edge  e  to  Edge  e',  and  Diffracted  From 
Edge  e'  to  the  Observation  Point. 


Observation 

Point 


Figure  13.  The  Triply-Reflected  Ray  Reflected  From  Plate  K  to  Plate  L,  Reflected  From  Plate  L 
to  Plate  M,  and  Reflected  From  Plate  M  to  the  Observation  Point. 
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Figure  14.  The  Reflected-Reflected-Diffracted  Ray  Reflected  From  Plate  K  to  Plate  L,  Reflected 
From  Plate  L  to  Edge  e,  and  Diffracted  From  Edge  e  to  the  Observation  Point. 


The  second  step  is  to  determine  whether  reflection  points  exist  within  Plates  K  and  L,  and  a 
diffraction  point  ( Q e)  on  Edge  c.  With  the  exception  of  the  case  of  an  edge  and  a  line  that  contains  the 
second  image  point  and  the  observation  being  parallel,  the  diffraction  point  always  exists  in  the  two- 
dimensional  case.  Qe  and  the  second  image  location  are  used  as  two  points  on  a  line.  The  reflection  point 
(RJL)  where  this  line  intercepts  Plate  L  is  found.  Then  using  RL  and  the  first  image  location  as  two  points 
on  a  line,  the  reflecting  point  RK,  which  is  the  interception  point  on  Plate  K,  is  calculated.  If  either  RL, 
RK,  or  Qc  do  not  exist,  the  RRD  ray  will  not  exist. 

The  third  step  is  to  test  whether  or  not  the  ray  from  the  source  to  RK,  RK  to  RL.,  RL  to  and 
Qe  to  the  observation  point  is  blocked  by  the  other  plates.  If  an  RRD  ray  passes  the  existence  tests,  its 
wave  amplitude  and  phase  are  calculated  and  added  to  the  previous  contributions.  This  algorithm  is 
performed  on  all  forward-traveling  ray  combinations  of  plates  and  edges.  The  final  result  becomes  the 
total  reflected-rcflected-diffracted  wave  field. 

Reflected-Diffracted-Reflected  Ray 

The  reficctcd-diffracted-reflected  (RDR)  ray  is  reflected  once  from  a  plate,  diffracted  by  an  edge 
that  is  not  part  of  the  plate,  and  reflected  once  more  from  another  plate.  Figure  15  depicts  an  RDR  ray. 

The  first  step  in  the  ray  existence  tests  is  to  determine  that  the  edge  (Edge  e  in  this  case)  is  located 
above  both  reflecting  plates  (Plates  K  and  M  are  considered).  If  this  condition  is  met,  the  source  and 
observation  point  are  imaged  with  respect  to  Plate  K  and  Plate  M,  respectively;  also,  an  image  of  Edge 
e  with  respect  to  Plate  K  is  calculated.  By  using  an  image  of  Edge  e,  the  image  source  and  the  image  of 
the  observation  point,  the  image  diffraction  point  (Qei)  is  determined.  The  method  in  determining  the 
diffraction  point  is  explained  previously.  Qei  is  then  imaged  again  with  respect  to  Plate  K  to  a  point  on 
Edge  e,  and  is  the  diffraction  point  (C^). 

The  second  step  is  to  test  for  the  existence  of  the  reflection  points  on  Plates  K  and  M.  The 
reflection  point  on  Plate  K  (RK)  is  determined  by  generating  a  line  that  contains  Qe  and  the  image  source, 
intercepting  Plate  K;  similarly,  the  reflection  point  on  Plate  M  (RM)  is  located  where  a  line  that  contains 
Qe  and  the  image  of  the  observation  point  intersects  Plate  M. 
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Figure  15.  The  Reflected-Diffracted-Reflected  Ray  From  Plate  K  to  Edge  e,  Diffracted  From  Edge 
e  to  Plate  M,  and  Reflected  From  Plate  M  to  the  Observation  Point. 
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Once  RK  and  RM  are  determined  to  lie  on  Plates  K  and  M,  the  last  step  required  is  the  blockage 
test.  The  line  of  sight  from  the  source  to  RK,  from  RK  to  from  Qe  to  RM,  and  from  RM  to  the 
observer  must  not  be  blocked  by  the  other  plates.  If  the  ray  passes  the  existence  test,  its  wave  amplitude 
and  phase  are  calculated  and  added  to  the  total.  This  process  is  repeated  until  all  forward  traveling  ray 
combinations  of  plates  and  edges  are  examined.  The  final  result  becomes  the  total  reflected-diffracted- 
reflected  field. 


Reflected-Diffracted-Dijfracted  Ray 

The  reflected-diffracted-diffracted  (RDD)  ray  is  defined  as  a  ray  reflected  from  a  plate  and  diffracted 
from  two  different  edges,  where  the  diffraction  edges  are  not  part  of  the  plate.  An  example  of  this  ray 
geometry  is  given  in  Figure  16. 

The  first  step  to  test  for  ray  existence  is  to  locate  the  image  source  with  respect  to  Plate  K.  Once 
this  image  point  is  found,  the  diffraction  points  on  Edges  e  and  e'  are  determined.  The  second  step  is  to 
determine  the  existence  of  the  reflection  point  (K  in  this  case).  This  is  done  by  generating  a  line  that 
contains  the  image  source  and  the  diffaction  point  (e  in  this  case)  and  testing  whether  it  intercepts  Plate 
K.  The  third  step  is  the  blockage  test.  A  line  that  originates  from  the  source  to  RK,  from  RK  to  Qe,  from 
Qe  to  Qe',  and  from  Qe'  to  the  observation  point  must  not  be  blocked  by  the  other  plates. 

Once  the  ray  passes  the  existence  tests,  its  wave  amplitude  and  phase  are  calculated  and  added  to 
the  previous  contributions.  This  algorithm  is  performed  on  all  forward-traveling  ray  combinations  of 
existing  plates  and  edges.  The  final  result  becomes  the  total  reflected-diffracted-diffractcd  field. 

Diffracted-Reflected-Reflected  Ray 

The  diffracted-refiected-reflected  (DRR)  ray  is  diffracted  once  from  an  edge  which  must  not  be  part 
of  a  reflecting  plate,  and  reflected  from  two  different  plates.  The  ray  geometry  is  illustrated  in  Figure  17. 

The  first  step  in  the  ray  existence  test  is  to  determine  that  the  diffraction  point  and  both  reflection 
points  exist.  This  is  done  by  reflecting  the  observer  location  with  respect  to  Plate  M.  and  further  reflecting 
this  image  with  respect  to  Plate  L  to  form  the  secondary-image  point.  From  the  secondary-image  observer 


Figure  16.  The  Reflected-Diffracted-Diffracted  Ray  Reflected  From  Plate  K  to  Edge  e.  Diffracted 
From  Edge  e  to  Edge  e',  and  Diffracted  From  Edge  e'  to  the  Observation  Point. 
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Figure  17.  The  Diffracted-Reflected-Reflected  Ray  Diffracted  From  Edge  e  to  Plate  L;  Reflected 
From  Plate  L  to  Plate  M,  and  Reflected  From  Plate  M  to  the  Observation  Point. 


point,  the  source,  and  the  Edge  e,  the  location  of  the  diffraction  point  (Qe)  is  calculated.  The  method  of 
obtaining  Qe  has  been  discussed  previously.  Qe  and  the  secondary-image  observer  location  form  a  line, 
whose  intersection  point  (RL)  on  the  extension  of  Plate  L  is  calculated;  similarly,  RL  and  the  first  image 
observation  point  form  another  line,  whose  intersection  point  (RM)  on  the  extension  of  Plate  M  is 
determined.  Both  RL  and  RM  must  lie  on  Plates  L  and  M  for  the  DRR  ray  to  exist. 

The  last  step  is  the  blockage  test.  Lines  from  the  source  to  Qc,  from  Qe  to  RL,  from  RL  to  RM, 
and  from  RM  to  the  observation  point  must  not  be  blocked  by  the  other  plates.  If  the  line  is  not  blocked, 
the  DRR  ray  wave  amplitude  and  phase  arc  calculated  and  added  to  the  previous  contributions. 

The  process  is  repeated  until  all  edge  and  plates  in  forward-direction  combinations  for  the  DRR  ray 
have  been  examined.  The  final  result  is  the  total  diffracted-reflected-reflected  field. 

Diffracted-Reflected-Diffracted  Ray 

The  diffracted-reflected -diffracted  (DRD)  ray  is  one  that  is  diffracted  from  an  edge,  reflected  from 
a  plate  that  is  not  part  of  the  same  edge,  and  diffracted  once  more  from  another  edge  which  is  not  part 
of  the  reflecting  plate.  An  example  of  this  ray  geometry  is  illustrated  in  Figure  18. 

The  first  condition  to  test  the  DRD  ray  is  possible  is  that  the  diffraction  points  on  both  edges  (e  and 
e'  in  this  case)  and  the  reflection  point  on  a  plate  (L  for  example)  must  exist.  One  method  is  to  reflect 
both  Edge  e'  and  the  observation  point  with  respect  to  Plate  L.  Then,  using  the  source.  Edge  e,  image 
of  Edge  e',  and  the  image  of  the  observation  point,  the  diffraction  points  are  determined  as  for  the  DD 
ray. 


The  second  step  is  to  test  whether  the  reflection  point  is  on  Plate  L.  The  reflection  point  is  obtained 
by  generating  a  line  that  contains  both  diffraction  points  and  intersects  Plate  I  at  RL.  For  the  DRD  ray 
to  exist,  RL  must  be  within  the  boundaries  of  Plate  L.  It  should  be  noted  that  one  diffraction  point  is  on 
the  image  of  Edge  e'.  Therefore,  before  starting  blockage  tests,  the  actual  diffraction  point  is  needed.  The 
resultant  point  (Q/)  is  thus  located  on  Edge  e'. 

The  final  condition  for  the  ray  to  exist  is  no  blockage.  The  lines  of  sight  from  the  source  to  Qe. 
from  Qt  to  RL,  from  RL  to  Q/,  and  from  Q/  to  the  observation  point  must  not  be  blocked  by  other  plates. 
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Figure  18.  The  Diffracted-Reflected-Diffracted-Ray  Diffracted  From  Edges  e,  to  Plate  L,  Reflected 
From  Plate  L  to  Edge  e',  and  Diffracted  From  Edge  e'  to  the  Observation  Point. 

If  not  blocked,  the  ray  wave  amplitude  and  phase  are  calculated  and  added  to  the  previous  contributions. 
The  process  is  repeated  until  all  edges  and  plates  in  forward-direction  combinations  for  the  DRD  ray  have 
been  tested.  The  final  result  is  the  total  diffracted-reflected-diffracted  field. 

Diffracted-Diffracted-Reflected  Ray 

The  diffracted-diffracted-reflected  (DDR)  ray  is  defined  as  a  ray  diffracted  from  two  different  edges 
and  reflected  once  from  a  plate  that  does  not  touch  those  diffracting  edges.  An  example  of  this  ray 
geometry  is  illustrated  in  Figure  19. 

The  first  step  is  to  locate  the  image  of  the  observation  point  with  respect  to  the  reflecting  plate 
(Plate  M  in  this  case).  The  image  location  must  lie  below  the  plane  that  contains  reflecting  Plate  M.  If 
this  condition  is  met,  the  next  step  is  performed. 

The  second  step  is  to  find  the  image  diffraction  points  on  both  edges  (Edges  e  and  e')  and  the 
reflection  point.  The  method  is  very  similar  to  the  one  used  for  the  doubly-diffracted  ray.  The  only 
difference  is  that  instead  of  the  real  observation  point,  the  image  observation  point  is  used.  Once  Qe  and 
Q/  are  found,  the  reflection  point  on  Plate  M  is  calculated.  Q/  and  the  image  observation  point  form  a 
line  that  intercepts  the  plane  that  contains  Plate  M.  If  the  reflection  point  (RM)  lies  within  Plate  M.  the 
blockage  test  is  performed. 

The  blockage  test  is  the  last  step  in  the  ray  existence  algorithm.  The  lines  of  sight  propagation  path 
from  the  source  to  Qt,  from  Qe  to  Qe\  from  Qe'  to  RM,  and  from  RM  to  the  observation  point  must  not 
be  blocked  by  the  other  plates. 

If  the  DDR  ray  passes  all  three  steps,  it  exists  and  the  wave  amplitude  and  phase  are  calculated  and 
added  to  the  previous  contributions.  This  algorithm  is  repeated  on  all  forward-direction  combinations  of 
plates  and  edges.  The  final  result  becomes  the  total  diffracted-diffracted-refiected  field. 

Triply-Diffracted  Ray 

The  triply-diffracted  (DDD)  ray  is  diffracted  by  three  different  edges.  An  example  of  this  ray  is 
illustrated  in  Figure  20.  Before  starting  the  blockage  tests,  the  three  diffraction  points,  Qe.  Q/.  and  Qe" 
are  located.  Then  blockage  tests  of  the  lines  of  sight  from  the  source  to  Qe,  from  Qe  to  Qe\  from  Q/  to 
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Figure  19.  The  Diffracted-Diffracted-Reflected  Ray  Diffracted  From  Edge  e  to  Edge  e',  Diffracted 
From  Edge  e"  to  Plate  M,  and  Reflected  From  Plate  M  to  the  Observation  Point. 
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Figure  20.  The  Diffracted-Diffracted-Diffracted-Ray  Diffracted  From  Edges  e,  e',  and  e"  to  the 
Observation  Point. 


Qe",  and  from  Qe"  to  the  observation  point  are  performed.  If  there  is  no  line  of  sight  blockage,  the  DDD 
ray  exists  and  its  amplitude  is  evaluated  and  added  to  the  sum. 


Description  of  Computer  Programs  for  the  Two-Dimensional  GTD  Model 

Most  of  the  two-dimensional  GTD  model  computer  subroutines  used  by  Luebbers  (1985)  are  valid 
and  used  in  this  three-dimensional  model.  In  this  section  we  give  a  brief  functional  explanation  of  the 
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most  important  subroutines.  In  the  next  four  sections  we  describe  the  revisions  and  improvements  of 
computer  programs  for  a  realistic  three-dimensional  model  of  SEL  calculation. 

Subroutine  VGROUN.f 

This  subroutine  VGROUN.f  initializes  all  of  the  values  required  for  the  GTD  calculations  and  acts 
as  a  buffer  between  the  main  program  and  the  peripheral  subroutines. 

The  observation  point  and  source  locations,  which  are  supplied  by  the  main  program,  are  organized 
for  efficient  manipulation.  Then  the  normal  vector  of  each  plate  is  determined  in  such  a  way  that  the  z- 
component  is  negative.  The  wedge  angle,  which  is  formed  by  two  connected  plates,  is  calculated  by  using 
their  normal  vectors  and  vector  cross-product  operation.  The  angle  is  measured  below  the  terrain  profile. 

To  reduce  the  calculation  time,  the  lengths  in  meters  are  converted  to  wavelengths.  After  con¬ 
version,  subroutines  VENGO.f  and  VETDIF.f  determine  the  relative  wave  field. 

Subroutine  VENGO.f 

The  subroutine  VENGO.f  calculates  the  total  wave  amplitude  and  phase,  which  result  from  the 
direct,  reflected  (R),  RR,  RD,  RRD,  RDR  and  RRR  rays.  The  subroutine  considers  the  direct  ray  first. 
The  direct  ray  originates  from  the  source  and  travels  directly  to  the  observation  point.  Subroutine 
VLOGIl.f  checks  any  possible  blockage.  If  there  is  no  blockage,  the  direct  wave  components  are 
calculated  by  subroutine  NWSRCE.f.  If  the  direct  ray  is  blocked,  the  direct  ray  contribution  is  zero. 

The  reflected  rays  are  considered  next.  AH  of  the  possible  singly-reflected  rays  from  every  plate 
are  examined  within  the  DO  statement  301,  where  LP  is  a  plate  index.  When  Plate  LP  is  under  test, 
subroutine  VEXIST  is  called  to  determine  the  image  source  (SI)  and  reflection  point  OQ  on  Plate  LP. 
If  OQ  does  exist,  the  blockage  algorithms  are  started  by  checking  for  the  blockage  from  ANT  (source 
point)  to  OQ,  and  from  OQ  to  FPTS  (observation  point).  If  no  blockage  occurs  (the  reflected  ray  exists), 
subroutines  NWSRCE  and  FCREFL  compute  the  reflected  wave  components  (ER)  from  Plate  LP.  Then 
ER  is  added  to  the  total  wave  (ESUM). 

After  the  reflected  ray  from  Plate  LP  is  evaluated,  the  doubly-reflected  (RR)  ray  is  considered.  The 
DO  statement  414  which  is  looped  inside  DO  statement  301,  will  examine  all  of  the  plates  LN  that  can 
possibly  interact  with  Plate  LP  to  contribute  a  RR  ray.  If  LN  is  equal  to  LP,  no  calculation  is  performed. 
Plate  LN  is  always  located  further  away  from  the  source  than  Plate  LP.  Therefore,  only  the  RR  rays 
traveling  in  the  forward  direction  are  included.  The  process  stops  after  the  last  plate  (NPLATE)  is 
considered. 

Subroutine  IMAGE  determines  a  secondary  image  location  (SlI)  of  the  image  source  (SI)  with 
respect  to  Plate  LN.  By  utilizing  pairs  (SII,  FPTS)  and  (SI,  ANT),  two  reflection  points  (OQ,  RP)  are 
determined  by  calling  VEXIST  twice.  OQ  and  RP  must  both  be  located  on  Plate  LN  and  LP,  respectively. 
If  they  are  not,  the  RR  ray  contribution  is  zero.  Also,  the  lines  from  ANT  to  RP,  RP  to  OQ,  and  OQ  to 
FPTS  must  not  be  blocked.  If  no  blockage  occurs,  routines  NWSRCE  and  FCREFL  are  called,  the  RR 
ray  components  from  Plate  LP  to  Plate  LN  are  determined  and  then  added  to  ESUM. 

During  the  RR  ray  determination,  several  parameters  required  for  the  RRD  ray  calculation  are 
calculated.  To  save  computer  time,  subroutine  VRRD  is  immediately  called  for  RRD  evaluations.  The 
output  from  subroutine  VRRD  is  added  to  the  total  wave  field  in  ESUM. 
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RD  ray  contributions  are  considered  next.  AH  possible  RD  ray  contributions  from  edges  are 
evaluated  with  DO  statement  426.  All  edges  past  the  edge  where  plates  LP+1  and  LP+2  are  connected 
are  considered.  The  RD  ray  is  reflected  from  Plate  LP  to  the  edge  and  diffracted  from  the  edge  to  the 
observation  point  (FPTS).  It  should  be  noted  that  DO  statement  426  is  nested  within  DO  statement  301. 

The  reflection  point  (RF1)  on  Plate  LP  and  the  possible  diffraction  point  (RP)  on  Edge  LM  are 
calculated  first.  If  both  RF1  and  RP  exist,  the  blockage  algorithms  are  performed  for  the  blockage  from 
ANT  to  RF1,  RF  to  RP,  and  RP  to  FPTS.  If  there  is  no  blockage,  routine  SDDIF  is  called  for  evaluation 
of  RD  ray  amplitude  determination.  The  wave  components  are  determined  by  routine  SDDIF  and  added 
to  the  total  field  in  ESUM. 

The  remaining  rays  considered  in  subroutine  VENGO  are  calculated  in  separate  subroutines  called 
by  VENGO.  Subroutine  VRDR  for  the  RDR  rays,  subroutine  VRDD  for  both  the  RDD  and  RDDD  rays, 
and  subroutine  VRRR  for  the  RRR  rays.  The  wave  components  for  each  ray  type  are  added  to  the  total 
in  array  ESUM. 

Subroutine  VETDIF.f 

Subroutine  VETDIF.f  calculates  seven  types  of  rays:  Singly-diffracted,  DR,  DD,  DDD,  DRD,  DDR 
and  DRR  rays.  It  is  called  only  from  subroutine  VGROUN. 

After  initialization,  the  edge  formed  by  two  connected  plates  (LP,  LP-1)  is  considered.  With  three 
inputs,  the  source  location,  the  edge,  and  the  observation  point,  subroutine  VBETA  is  called  for 
determining  the  diffraction  point  (DPT)  on  the  edge.  Then,  using  ANT,  DPT,  and  FPTS,  subroutine 
SDDIF  calculates  wave  components  of  the  singly-diffracted  ray.  The  components  are  then  added  to  the 
previous  contributions. 

The  DR  ray  is  considered  next,  and  it  must  be  located  on  plate  LJ.  Reflection  point  C  is  determined 
by  calling  subroutine  VEXIST.  If  either  DPT  or  C  do  not  exist,  the  DR  ray  contributions  are  set  to  zero, 
and  plate  LJ+1  is  considered.  Otherwise,  the  blockage  algorithm  begins  by  testing  the  lines  from  ANT 
to  DPT,  from  DPT  to  C,  and  from  C  to  FPTS.  If  there  is  no  blockage,  the  DR  ray  exists,  and  subroutine 
SDDIF  determines  the  ray  component. 

Following  DR  ray  calculations,  DRR  rays  are  examined.  This  is  done  by  simply  calling  subroutine 
VDRR.  The  wave  amplitudes  and  phases  of  DRR  rays  are  added  to  the  total. 

The  next  ray  to  be  considered  is  the  DRD  ray.  If  a  DRD  ray  exists,  the  edge  formed  by  connected 
plates  LR  and  LR-1,  where  Plate  LR  is  different  from  plates  LP  and  LJ,  is  considered.  It  should  be  noted 
that  LR-1  is  always  greater  than  LJ.  The  edge  formed  by  plates  LR  and  LR-1  and  FPTS  are  both  imaged 
with  respect  to  extension  of  Plate  LJ.  Using  the  source  location,  the  edge  (formed  by  plates  LP  and  LP-1), 
the  image  of  the  edge  (formed  by  plates  LR  and  LR-1),  and  the  image  of  FPTS,  the  diffraction  points  DPI 
and  DP2  are  calculated  by  subroutine  V2DIF.  If  DPI  and  DP2  exist,  subroutine  VEXIST  calculates  the 
reflection  point  location  DPT  on  Plate  LJ.  If  DPT  does  not  lie  on  Plate  LJ.  the  next  edge  is  considered. 
Otherwise,  the  blockage  algorithm  tests  for  blockage  of  the  lines  from  ANT  to  DPI,  from  DPI  to  DPT, 
from  DPT  to  DP2,  and  from  DP2  to  FPTS.  If  there  is  no  blockage,  subroutine  SDDIF  is  called  twice  for 
DRD  ray  calculations.  The  ray  contributions  are  added  to  the  total. 

The  next  ray  considered  in  subroutine  VETDIF  is  the  DD  ray.  This  ray  is  diffracted  by  two  edges. 
The  first  edge  is  formed  by  plates  LP  and  LP-1,  whereas  the  second  edge  is  formed  by  plates  LM  and 
LM-1.  Subroutine  V2DIF  determines  two  diffraction  points  (DPI,  DP2).  The  blockage  algorithm  tests 
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the  lines  from  ANT  to  DPI,  from  DPI  to  DP2,  and  from  DP2  to  FPTS.  If  there  is  no  blockage, 
subroutine  SDDIF  determines  DD  ray  amplitudes  and  phases,  and  they  are  added  to  the  total. 

The  next  ray  considered  is  the  DDR  ray.  The  reflecting  plate  is  LS.  The  process  begins  as  FPTS 
is  imaged  with  respect  to  Plate  LS.  It  is  necessary  that  the  image  must  be  located  below  the  plane  of  Plate 
LS.  Subroutine  VEXIST  is  called  to  calculate  the  reflection  point  C  on  Plate  LS.  If  C  does  not  exist,  the 
next  Plate  LS+1  is  considered.  Otherwise,  the  blockage  is  tested  on  the  lines  from  ANT  to  DPI,  from 
DPI  to  DP2,  from  DP2  to  C,  and  from  C  to  FPTS.  If  the  lines  are  not  blocked,  subroutine  SDDIF  is 
called  twice  for  DDR  ray  amplitude  and  phase  calculations,  and  the  results  are  added  to  the  corresponding 
total  field  components. 

Finally,  the  DDD  ray  is  considered.  Blockage  and  existence  tests  already  performed  for  the  DD  ray 
are  used.  Since  the  ray  can  propagate  from  edges  following  plates  LP  and  LM,  existence  and  blockage 
tests  are  done  from  Edge  LM  to  LS  and  LS  to  FPTS.  If  the  DDD  ray  can  exist,  it  is  evaluated  by  calling 
SDDIF  three  times  and  adding  the  result  to  the  total. 

Subroutines  VAMPHA.f  and  VAMPH3.f 

These  subroutines  convert  real  and  imaginary  components  of  complex  wave  fields  into  amplitudes 
and  phases. 

Subroutine  VMERGH.f 

This  subroutine  tests  for  ray  blockage  when  the  end  points  are  situated  on  two  consecutive  plates, 
K  and  L.  It  is  assumed  that  Plate  L  is  located  further  away  from  the  source  than  Plate  K. 

Subroutine  VECTOR.f 

This  can  perform  six  vector  operations,  which  are:  (1)  cross  product;  (2)  normalized  cross  product; 
(3)  unit  vector  of  a  difference  between  two  vectors;  (4)  dot  product;  (5)  projection  of  two-directional 
vectors;  (6)  vector  difference  between  two  position  vectors. 

More  detailed  descriptions  about  the  subroutines  of  the  two-dimensional  GTD  model  can  be  found 
in  Luebbers  (1985).  Documentation  and  help  are  also  included  in  the  computer  programs. 
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3  THE  PROPAGATION  MODEL 


From  Two  Dimensions  to  Three  Dimensions 

Here  the  geometrical  theory  assumption  to  the  reflection  of  a  sound  wave  by  a  plate  is  applied  using 
the  same  theory  for  either  two-  or  three-dimensional  terrain.  However,  the  diffraction  from  a  hill  of  finite 
length  is  most  different  from  the  diffraction  of  a  hill  of  infinite  length  in  cases  where  lower-frequency 
sound  waves  are  used.  Although  some  investigations  that  extend  GTD  to  lower  frequencies  could  have 
been  used,  it  was  expedient  for  the  present  study  to  continue  with  the  same  reflection  and  diffraction 
formulae.  In  this  regard,  the  GTD  method  itself  is  somewhat  crude,  but  one  of  the  best  models  for  this 
situation  short  of  a  brute-force  numerical  procedure. 

The  amount  of  diffracted  energy  re-radiated  from  an  edge  is  determined  by  the  diffraction  coeffi¬ 
cient,  which  is  a  function  of  the  wedge  angle  and  the  incident  and  diffracted  ray  geometries.  The  detailed 
deduction  of  formulas  used  for  two-dimensional  computations  was  given  by  Luebbers  for  single  wedge 
diffraction.  The  formulas  for  double  wedge  diffraction  are  given  in  Chapter  4. 

Subroutine  VGROUN  in  the  original  computer  program  acts  as  a  buffer  between  the  main  GRASS 
program  and  the  peripheral  subroutines.  In  VGROUN,  the  x-component  of  the  plates  (where  the  x-axis 
is  parallel  to  the  hills)  is  extended  from  -50,000  to  50,000  ft.  Normal  vectors  and  wedge  angles  of  the 
plates  are  calculated  from  the  four  comer  positions  of  each  plate.  A  three-dimensional  algorithm  was 
developed  to  determine  the  normal  vectors  and  wedge  angles  of  the  plates  to  fit  any  arbitrarily  oriented 
connected  plates.  It  is  therefore  more  convenient  to  directly  specify  the  position  of  the  comers  as  input 
data  to  represent  the  terrain. 

The  actual  calculation  of  the  sound  pressure  is  conducted  in  the  subroutines  NWSRCE  and  SDDIF. 
Although  subroutine  VGROUN  returns  the  sound  pressure,  the  sound  pressure  for  each  ray  is  reported  by 
subroutine  VAMPH3.  Other  related  subroutines  include  VENGO,  VENDIFF,  VRDD,  VRDR,  VRRD, 
VMERGH,  and  VLOGI1.  The  subroutines  were  converted  or  modified  for  predictions  of  the  sound  field 
in  three  dimensions  above  arbitrary  three-dimensional  terrain.  The  more  complex  changes  are  explained 
in  the  following  chapters. 


Frequency  Switch  and  Blast  Waves 

The  GTD  propagation  model  was  used  here  to  calculate  the  noise  levels  produced  by  heavy  weapons 
that  generate  blast  waves.  Four  parameters  are  usually  used  to  characterize  a  blast  wave:  peak  pressure, 
duration,  central  frequency,  and  half-width  of  the  frequency  spectrum.  For  blasts  produced  by  Army 
weapons,  it  was  necessary  to  consider  a  total  frequency  band  from  one-third  octave  band  10  (10  Hz)  to 
one-third  octave  band  20  (100  Hz).  For  each  one-third  octave  band,  the  propagation  effects  on  only  the 
band  center  frequency  were  considered.  In  this  way,  a  wave  energy  evolution  very  near  that  of  real, 
measured  waves  can  be  obtained  without  performing  calculations  at  hundreds  of  frequencies. 

The  user  may  set  a  switch  index  (K,)  to  choose  one  of  five  options  for  the  frequency  characteristics 
of  the  sound  source.  When  K,  =  1,  a  single  frequency  source  is  selected.  For  single  frequency 
calculations,  the  user  inputs  the  frequency  and  the  pressure  at  unit  distance  (1  m)  from  the  source.  When 
Kf  =  2,  multiple  frequencies  arc  used.  The  user  inputs  all  of  these  frequencies  and  the  pressure  values  of 
the  respective  frequencies  for  a  specified  unit  distance  from  the  source.  Kf  =  3  is  selected  for  blast  wave 
calculations.  For  blast  calculations,  the  user  inputs:  a  maximum  pressure  (Pm„);  the  duration  (Td);  the 
number  of  frequencies  (N,rc),  the  central  frequency,  and  half-width,  in  response  to  the  prompt  by  the 
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program.  Kf  =  4  is  also  used  for  blast  waves  but  in  this  case  the  user  need  only  input  charge  and 
frequency  parameters.  The  computer  will  calculate  the  other  blast  wave  parameters  Pmax  and  Td  from  the 
charge  using  the  following  default  algorithm: 


Pmx  -  7-44X104 


f  w  ^°-4 
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[Eq  3] 
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lEq  4] 


where  the  charge  weight  W  is  in  pounds,  Pmax  is  the  maximum  pressure  in  pascals,  and  Td  is  the  duration 
in  seconds.  Kf  =  5  is  also  used  for  blast  waves,  but  for  this  index  the  only  input  is  the  charge  size  of  the 
weapon.  For  Kf  =  5  the  computer  will  not  only  calculate  ,x  ar  „  Td  but  also  set  the  frequency  range 
to  vary  from  one-third  octave  band  10  to  one-third  octa  -and  20.  Subroutines  blast.f  and  blastl.f 
calculate  the  sound  pressures  at  each  frequency  one  unit  distance  from  the  explosion. 


Noise-Level  Weighting  Options 

Because  the  sensitivity  of  the  human  ear  varies  with  frequency,  the  significance  of  noise  levels  at 
different  frequencies  also  varies.  The  International  Standard  for  Sound  Level  Meters,  DEC  651-1979,  and 
its  revision,  the  American  National  Standard  Institute  (ANSI)  Specification  for  Sound  Level  Meters ,  ANSI 
SI. 4-1983,  provide  a  set  of  A-,  B-,  C-weighting  networks  for  sound-level  meters.  As  a  reference,  the 
revised  document  also  gave  similar  expressions  for  the  D-  and  E-  weighting  networks. 

A  switch  is  set  in  our  computer  program  for  the  choice  of  frequency  weighting.  For  the  choice 
Kw  =  0,  no  weighting  is  used.  Kw  =  1,  2,  3, 4,  and  5  use  A-,  B-,  C-,  D-,  and  E- weightings,  respectively. 
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4  OTHER  PROGRAM  IMPROVEMENTS 


Subroutine  V2DIF  determines  two  diffraction  points  that  exist  on  two  different  edges.  The 
arguments  included  in  the  original  codes  included  only  one  diffraction  angle.  This  works  when  the  two 
edges  are  parallel,  but  for  three-dimensional  terrain,  the  two  edges  may  not  be  parallel.  Two  diffraction 
angles  were  set  up  for  the  diffractions  from  the  two  edges  (Figure  21).  In  the  original  codes,  a  plane  that 
includes  two  edges  is  introduced  as  the  reference  plane.  In  three  dimensions,  the  two  wedges  may  not 
necessarily  be  on  one  plane.  Figure  21  shows  A  as  the  source  position  and  F  as  the  observer  position. 
The  unit  vector  u,  is  along  the  first  wedge  £,  £2;  the  unit  vector  is  along  the  second  wedge  £3  -  £4; 

(?,  and  Q2  are  two  assumed  diffraction  points.  Line  AO  is  perpendicular  to  u,  in  the  plane  composed  by 
A  and  vector  Line  £)2F'  is  perpendicular  to  m,  in  the  plane  composed  by  Q2  and  Using  single 
diffraction  geometry,  the  relation  between  r,,  r2,  r3,  and  x,  is: 


A  similar  relation  exists  between  r,,  r2,  r3,  and  x2  for  a  single  diffraction  of  the  second  wedge. 

Therefore,  given  a  source  position  A  and  one  diffraction  point  Q2,  another  diffraction  point  Q,  can 
be  calculated:  and  given  a  diffraction  point  Q,  and  an  observation  position  F,  the  other  diffraction  point 
Q2  can  be  obtained.  Therefore,  given  A  and  F,  there  are  two  equations  for  Q,  and  Q2,  making  the 
calculation  of  Q,  and  Q2  feasible.  The  new  subroutine  V2DIF  uses  an  iteration  method  to  solve  for  Q, 
and  Q2.  Test  runs  show  the  iteration  quickly  converges  with  a  criteria  of  one-percent  accuracy. 


Figure  21.  Geometry  of  Double  Diffraction. 
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The  subroutine  VMERGH  is  called  when  ray  blockage  occurs  where  ray  end  points  are  situated  on 
two  consecutive  plates,  1  and  2,  and  where  normal  vectors  N,  and  N2  point  upward  as  positive.  An 
alternative  method  now  fulfills  this  function  (Figure  22).  If  the  ray  intersects  the  plates  at  P,  and  P2  and 
the  dot  product  of  the  vector  P,  to  P2  with  N,  or  the  vector  P2  to  P,  with  N2  is  less  than  zero,  the  ray  will 
not  be  blocked  by  the  valley  edge  of  the  two  plates.  If  the  dot  product  is  larger  than  zero,  the  ray  will 
be  blocked  by  the  peak  edge  of  the  two  plates.  This  criterion  is  simpler  than  the  original  one, 
(Houshmand  1990)  and  the  user  need  not  assume  that  Plate  2  is  located  farther  from  the  source  than  Plate 
1,  as  in  the  original  subroutine. 

Two  bridges  were  created  between  “gtd”  and  "GRASS":  a  subroutine  “tunder.f’  was  written  in 
Fortran  language;  and  the  equivalent  function  “input.c"  was  written  in  C  language,  and  is  included  as  part 
of  this  combined  program  system.  Subroutine  “tunder.f’  is  used  to  find  the  plates  that  are  under  the  line 
connecting  the  source  point  and  observation  point.  This  subroutine  also  prepares  formatted  terrain  data 
on  the  positions  of  source,  the  observation  point,  and  the  four  comers  for  each  plates. 

For  accuracy,  the  subroutine  “tunder.f'  first  distinguishes  four  different  situations  according  to  the 
observation  position  with  respect  to  the  source  point.  If  the  angle  between  the  Source-Observer  Line 
(SOL)  and  positive  x-axis  is  less  than  or  equal  to  45  degrees,  the  geometrical  situation  is  called  a 
quadrant- 1  position.  If  the  angle  between  the  SOL  and  positive  y-axis  is  less  than  45  degrees,  the  observer 
is  said  to  be  in  quadrant  2  with  respect  to  the  source.  If  the  angle  between  the  SOL  and  the  negative 
x-axis  is  less  than  or  equal  to  45  degrees,  it  is  said  to  be  a  third-quadrant  case.  If  the  angle  between  the 
SOL  and  the  negative  y-axis  is  less  than  45  degrees,  the  observer  is  said  to  be  in  the  fourth  quadrant. 

The  variable  i  is  used  as  the  row  index  of  the  terrain  grid  cells,  and  j  as  the  column  index  (Figure 
23).  At  each  grid  point,  (ij),  GRASS  gives  the  elevation  value.  Let  the  bottom-left  comer  of  the  terrain 
be  the  Cartesian  origin  (x,,.  y0).  The  largest  row  index  is  “iT.”  The  largest  column  index  is  “/T.”  Figure 
23  gives  a  quadrant  1  situation.  Point  P  is  a  given  source  position  in  the  cell  with  a  lower-left  comer  of 
A  (js,  y'J.  The  position  (is,  jj  is  referred  to  as  the  source  comer.  Point  F  is  a  given  observation  position, 
which  is  in  the  cell  with  a  lower-left  comer  of  B  (iD,  jQ).  We  call  position  (r0,  ja)  the  observer  comer.  Line 
AB  is  the  SOL.  If  the  intersection  point  of  the  SOL  and  vertical  grid  line  of  column  index  yc  is  Q.  then 
the  nearest  lower  comer,  G,  is  called  the  crossing  comer  of  yc.  The  row  index  of  the  crossing  comer  is 
ic,  making  the  crossing  comer  of  the  source  the  source  comer,  and  the  crossing  comer  of  the  observer  the 


Figure  22.  Geometry  of  Testing  Blockage  of  Ray  P1P2  by  Wedge. 
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observer  comer.  Next  determine  the  bottom  comer  and  top  comer.  All  the  column  grid  lines  which  cross 
with  the  SOL  will  be  the  edges  of  plates.  According  to  our  convention,  the  number  of  edges  are  counted 
from  the  source  to  observer  comer.  Next  the  bottom  comer  (ib,  jb)  and  too  comer  (/,,  jt)  of  the  edges  are 
determined.  If  the  ic  of  the  farther  column  grid  is  larger  than  the  column  under  consideration,  then  ib 
should  be  set  equal  to  ic  +  1  and  the  top  comer  should  be  set  equal  to  rc  -  1.  If  the  ic  of  the  farther 
column  grid  is  the  same  as  the  present  column  grid,  then  set  i'b  =  ic  and  i,  =  ic  -  1.  If  the  farther  column 
has  a  smaller  /c,  the  ib  should  equal  ic  but  set  it  =  ic  -  2.  Another  principle  to  determine  /„  and  /,  is:  if  the 
SOL  crosses  precisely  at  /b,  increment  /b  by  1,  i.e.,  move  down  one  row;  if  the  SOL  crosses  precisely  at 
i„  then  decrease  i,  by  1,  i.e.,  move  up  one  row.  This  is  necessary  to  avoid  shifting  the  reflection  point  or 
diffraction  point  outside  the  plate  or  wedge.  The  last  principle  prevents  any  i  from  being  larger  than  /T 
or  less  than  I,  and  similarly  for  j.  According  to  our  convention,  the  comers  should  be  ordered  from 
bottom  comer  to  top  comer  for  the  odd  edges  but  from  the  top  to  the  bottom  for  the  even  edges,  as 
demonstrated  by  the  arrows  in  Figure  23. 

Comers  c„  c2,  c3,  and  c4  compose  the  first  plate.  Comers  c3,  c4,  c5,  and  c6  compose  the  second  plate 
and  so  on.  Because  the  elevations  of  the  comers  are  from  terrain  data  generated  by  GRASS,  it  is  not 
guaranteed  that  the  four  comers  of  each  plate  will  lie  in  one  plane.  But  within  the  resolution  defined  here, 
it  is  reasonable  to  force  each  plate  to  fit  in  one  plane.  The  computer  “adjusts”  the  elevation  of  c4  to  make 


Figure  23.  Method  of  Finding  and  Rearranging  Plates  Under  the  Source-Observer  Line. 
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c4  into  the  last  comer  of  the  plane  of  the  first  plate  composed  by  the  comers  c,.  c2,  and  c3.  The  computer 
“adjusts”  the  elevation  of  the  last  comer  of  the  second  plate,  c6,  to  lie  in  the  plane  formed  by  the  adjusted 
new  comer  c4  and  the  original  comers  c3  and  c5  and  to  assure  that  the  neighboring  plates  have  a  common 
wedge.  In  this  way  the  elevation  of  the  bottom  comers  of  the  even  wedges  and  the  top  comers  of  the  odd 
wedges  will  be  adjusted. 

By  the  same  reasoning,  the  computer  “adjusts”  the  elevation  of  the  source  position  to  ensure  that 
the  source  is  situated  2.5  m  above  the  source  column  (for  quadrants  I  and  3)  or  source  row  (  for  quadrants 
2  and  4  ).  The  computer  also  “adjusts"  the  observer's  elevation  to  ensure  that  the  observer  is  positioned 
1.25  m  above  the  observer  column  (  for  quadrants  1  and  3  )  or  observer  row  (  for  quadrants  2  and  4  ). 
The  calculations  and  coordinates  are  different,  but  the  principles  arc  the  same  in  all  four  quadrants. 

The  GTD  program  uses  conventional  positions  to  determine  x  and  y  coordinate  axes.  In  quadrant 
1.  (-y)  of  the  terrain  axis  is  used  as  (+x)  of  the  “gtd"  program  axis,  and  (+x)  of  the  terrain  axis  as  (+y) 
of  the  “gtd”  axis.  In  quadrant  2,  (+x)  of  the  terrain  axis  is  also  taken  as  (+x)  of  the  “gtd”  program  axis, 
and  the  (+y)  of  the  terrain  axis  is  also  taken  as  (+y)  of  the  “gtd”  axis.  In  quadrant  3,  the  (+y)  direction 
of  the  terrain  is  the  (+x)  of  “gtd"  program  direction,  and  the  (-x)  of  terrain  will  be  the  (+y)  of  the  “gtd” 
program  axis.  In  quadrant  4,  the  (-y)  of  terrain  will  be  the  (+y)  of  the  “gtd"  axis,  and  the  (-x)  of  terrain 
will  be  the  (+x)  of  the  “gtd"  axis. 

For  the  sake  of  efficiency,  the  first  step  is  to  check  neighboring  plates  to  see  if  they  form  very  small 
edge  angles.  If  so,  they  are  replaced  with  larger  plates  before  re-arranging  the  comer  data.  For  this  study, 
an  angle  of  less  than  0. 1  radian  was  set  as  the  small-angle  criterion.  All  of  the  documented  subroutines 
are  listed  in  the  Appendix  to  this  report. 
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5  THE  OVERALL  CALCULATION 


As  mentioned  in  Chapter  1 ,  terrain  data  for  the  calculation  of  SEL  can  either  be  given  by  the  user 
or  input  from  GRASS  output.  When  the  data  arc  given  by  the  user,  the  computer  program  will  ask  the 
user  to  input  a  file  name  which  contains  terrain  data  and  then  conduct  the  calculation.  Here,  conveniently 
terrain  data  arc  supplied  by  an  existing  GRASS  data  base. 

To  obtain  an  SEL  map  as  simply  as  possible,  all  commands  are  collected  in  a  profile  called 
"gtdcom,”  which  is  saved  in  (USACERL)  directory  cer!:/en/usr/zhuang/gtd  (Figure  24).  Note  that 
directory  names  are  given  here  only  for  an  example.  The  user  may  modify  “gtdcom”  to  specify  any 
directories  that  contain  the  programs.  A  user  in  a  terminal  connected  to  the  CERL  computer  need  only 
type  the  command  “gtdcom”  along  with  the  following  input  parameters: 

1.  x,y  coordinates  of  the  source  position 

2.  Weighting  index 

3.  Charge  size  in  pounds. 

Before  running  “gtdcom,”  the  user  must  specify  the  exact  terrain  data  parameters.  The  terrain  data 
should  already  be  located  in  a  GRASS  data  base.  The  geographic  region  for  the  terrain  data  should  be 
correct.  The  resolution  of  the  cells  should  be  appropriate,  i.e.,  not  so  fine  that  the  number  of  cells  is  too 
large  (usually  a  region  no  larger  than  50  by  50)  and  also  not  so  rough  that  they  provide  poor  accuracy. 
All  these  items  can  be  checked  by  logging  in  to  any  GRASS  computer,  and  by  using  related  tools 
provided  in  GRASS  (Wcstcrvelt  1989). 

After  the  computer  receives  these  input  parameters,  it  will  follow  the  procedure  shown  in  Figure  24. 
The  USACERL  computer  first  connects  to  the  computer  (zorro)  on  which  GRASS  terrain  data  and 
programs  are  stored.  This  computer  will  enter  the  user’s  GRASS  environment,  write  the  geographic  region 
definition  parameters  to  the  terrain  file  "terra.in”  as  the  header,  write  the  raster  (cell)  elevation  data  to  the 
terrain  file,  and  remotely  copy  the  terrain  file  back  to  the  CERL  computer  and  store  it  in  the  directory 
cerl:/en/usr/zhuang/gtd/cases/tcrra.in  for  the  calculation  of  “gtd.”  The  CERL  computer  will  disconnect 
from  the  other  machine  and  run  “gtd.”  (The  group  of  programs  comprising  “gtd”  is  shown  in  the  figure 
inside  the  square  solid  line.)  The  gtd  program  will  usually  take  several  hours  to  finish  running,  depending 
on  the  number  of  cells.  One  value  of  the  SEL  for  each  cell  is  calculated  by  "gtd”  computer  codes.  It 
takes  about  10  seconds  for  the  present  CERL  computer  to  compute  a  single  value.  The  CERL  computer 
will  store  the  SEL  data  in  a  file  named  “final.dat.” 

After  obtaining  all  the  SEL  values,  the  USACERL  computer  will  remotely  copy  the  file  “final.dat” 
to  zorro  and  reconnect  to  the  zorro  computer.  Then  the  GRASS  tool  “r.in.ascii”  will  be  used  to  convert 
ASCII  “final.dat”  file  into  a  binary  raster  file  "final.map.”  Next,  a  GRASS  tool  “p.map”  will  be  used  to 
draw  a  hardcopy  SEL  map. 

The  three  main  processes  of  the  calculation  are  to:  (1)  get  terrain  data  from  GRASS,  (2)  run  “gtd” 
to  get  SEL  values,  and  (3)  display  the  contour  map  of  SEL.  Three  by-product  files  “sumray.dat,” 
“eachry.dat”  and  "tcst.dat"  may  be  generated  during  the  calculation  of  “gtd.”  The  file  "sumray.dat" 
contains  more  data  of  each  SEL  calculation,  including  the  summation  of  all  possible  rays  (direct  ray, 
reflected  rays,  and  diffractcr  rays)  without  weighting  for  each  frequency.  The  file  "eachry.dat”  contains 
the  calculated  results  of  sound  amplitude  without  weighting  for  each  ray  of  sound  waves.  The  file 
“test.dat”  contains  necessary  debugging  data.  If  the  user  is  interested  in  obtaining  these  files,  the 
comments  must  be  cancelled.  The  two  files  can  be  used  in  checking  or  examining  the  numbers  in 
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Figure  24.  Program  Flow  Chart. 


“final.dat."  In  the  appendix  program  listing,  the  subroutine  names  with  extention  “.w”  are  for  writing  out 
all  results  that  fill  the  files  “sumray.dat,”  “eachry.dat,”  and  “test.dat”  to  standard  output. 

There  are  two  interfaces  between  the  “gtd”  and  “GRASS”  systems:  the  interface  with  the  terrain 
data,  and  the  routine  that  draws  the  SEL  map.  These  interfaces  correspond  with  two  corresponding 
commands  on  the  GRASS  computer:  “tcoman”  and  “tdisp.” 

These  computer  programs  implicitly  possess  the  ability  to  deal  with  more  flexible  situations.  In 
addition  to  the  three  switches  inside  the  program  “maiaf  ’  (frequency  switch,  weighting  switch,  and  ray 
control  index),  the  user  can  input  arbitrary  terrain  data  by  creating  the  file  “terrain.”  To  do  this,  the  user 
must  delete  the  commands  in  the  command  file  “gtdcom”  that  generate  a  “terra.in”  file. 

Even  though  the  “gtd”  program  sets  the  maximum  number  of  frequencies  to  50,  this  study  used  only 
frequencies  from  10  Hz  to  100  Hz,  to  save  computer  time,  and  to  provide  enough  range  to  approximate 
heavy  weapon  blast  waves.  The  user  may  easily  change  these  numbers  to  fit  different  puiposes. 

Pressure  is  calculated  as  the  summation  of  the  following  rays:  (1)  direct  ray,  (2)  reflected  ray,  (3) 
doubly-reflected  ray,  RR,  (4)  RD  (reflected  then  diffracted  ray),  (5)  RRD,  (6)  RDR,  (7)  D  (diffracted  ray), 
(8)  DR,  (9)  DRD,  (10)  DD,  (11)  DDR,  (12)  DDR,  (13)  RDD,  (14)  RRR,  (15)  DDD,  (16)  RDDD.  These 
16  numbers  determine  whether  to  include  or  exclude  the  relevant  ray.  If  a  number’s  value  is  zero,  the 
relevant  ray  will  be  included,  but  if  a  number’s  value  is  1,  the  relevant  ray  will  be  excluded.  Because  the 
triple  diffraction  ray  is  usually  very  weak,  and  because  it  was  necessary  to  save  computer  time,  rays  15 
and  16,  which  have  DDD  rays,  were  simply  dropped.  In  the  final  version  of  programs  for  this  study,  all 
other  14  ray  indexes  were  set  to  zero. 
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6  TWO  SAMPLE  RUNS 


The  GRASS  data  file  “Pinyon"  shows  typical  terrain  features  at  Pinyon  Canyon,  CO.  Figure  25 
shows  sample  elevation  data,  smoothed  from  the  original  50-meter  resolution  geographical  data.  Before 
running  the  programs,  the  terrain  data  need  some  preparation: 

1.  Enter  the  GRASS  environment  and  choose  “mapset  pinyon” 

2.  Use  the  GRASS  tool  “r.neighbors”  to  smooth  the  elevation  data 

3.  Use  the  “g. region  tool  to  set  resolution  =  4000  (meters),  and  the  east-west-north-south 
boundaries.  The  job  may  now  be  run,  as  explained  in  the  previous  section.  Figure  25  shows  the  contents 
of  the  intermediate  file  “terra.in.”  The  format  of  file  ”terra.in”  is  specified  by  “gtd”  programs.  The  first 
four  lines  give  the  boundary  coordinates  of  the  terrain  area  considered.  The  fifth  line  gives  the  total  row 
number  “/T.”  The  sixth  line  is  the  total  column  number  “yT."  The  matrix  takes  the  form  rT  by  yT,  which 
gives  the  elevation  of  each  pixel  in  the  area.  Therefore,  the  first  line  of  the  matrix  is  for  row  1,  column 
1,2,  3...  yT. 

Figure  26  shows  the  contour  map  of  elevations  of  the  terrain  in  the  area  described  in  Figure  25. 
The  upper  left  comer  is  the  highest  region,  and  the  bottom  right  comer  is  the  second  highest  hill.  The 
lowest  place  is  located  in  the  upper  right  comer.  Figure  27  shows  the  contour  map  of  a  calculated  SEL, 
the  final  result  generated  by  “gtd”  and  GRASS.  The  maps  shown  in  Figures  26  and  27  use  the  same 
region  and  grid  size.  The  SEL  mapped  in  Figure  27  is  from  a  1-lb  (0.453  kg)  explosive  source  at  the 
center  cell.  The  “final.dat”  file  contents  are  shown  in  Figure  28. 

SEL  field  measurements  were  taken  at  the  Fort  Leonard  Wood,  MO  training  site  in  1989.  Elevation 
data  are  not  available  for  this  region  in  a  GRASS  data  base.  Table  2  gives  estimated  elevation  data  taken 
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Figure  25.  Sample  Grid  Cell  Elevation  Values  at  Pinyon  Canyon,  CO 


39 


window:  north:  4159950  m 
south:  4127950  m 
west:  569950  m 

east:  699950  m 

grid  size:  40km  x  40km 


Figure  26.  Elevation  Contours  in  the  Pinyon  Canyon  Region. 
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Figure  27.  Contour  Map  of  Calculated  SEL  Values  for  Pinyon  Canyon,  CO. 


by  hand  from  a  contour  map.  Table  3  lists  the  positions  of  two  blast  sources  and  seven  microphones. 
Tables  4  and  5  give  calculated  SEL  values  for  the  training  site,  for  a  grid  size  of  1000  m.  Figures  29  and 
30  respectively  show  elevation  contours  based  on  Table  2  and  the  calculated  SEL  based  on  Tables  4  and 
5.  Both  Figures  29  and  30  are  based  on  the  following  parameters:  (1)  region:  north  78500m;  south 
72500m;  east  86500m;  west  68500m,  (2)  grid  size:  100m  x  100m.  Figure  29  was  derived  from 
information  taken  from  blast  1,  and  Figure  30  from  blast  2.  Table  6  lists  observed  SEL  data,  averaged 
from  original  measured  data.  The  numbers  in  the  parentheses  arc  root  mean  square  errors.  The  calculated 
values  are  the  SEL  values  multiplied  by  10,  to  create  three-digit  measurements,  and  to  meet  the  GRASS 
requirement  that  data  for  map  drawing  be  integers.  A  comparison  between  Tables  4  and  5  with  6  and  7 
shows  some  agreement  between  calculated  and  observed  SEL  variation  and  distance.  In  general,  the 
calculated  values  are  larger  than  the  measurements.  This  is  probably  because  the  dissipation  of  the 
atmosphere  and  the  roughness  of  the  ground  were  not  included  in  the  model.  In  addition,  other  errors  may 
have  arisen  from  the  elevation  estimation  and  the  large  grid  size. 

We  have  chosen  the  above  two  sites  because  they  represent,  in  our  opinion,  extreme  examples  of 
the  kind  of  terrain  that  may  be  encountered  at  Army  bases  around  the  world.  Fort  Leonard  Wood  site  has 
very  large  variations  in  altitude  from  place  to  place.  It  is  typical  of  the  terrain  in  Germany  where  there 
are  highly  visible  noise  problems. 
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Figure  28.  Sample  Grid  Cell  Values  (xlO)  for  Calculated  SELs  in  the  Pinyon  Canyon  Region. 
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Table  2 


Estimated  Elevation  Data  for  Fort  Leonard  Wood  Region 


(From  sheet  7559  II  Series  v  779  Ed  6-DMA  taken  from  10  m  contour 
lines  on  kilometer  grid  scale  1:50000  ) 


N,  E 

N,  E 

N, 
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N, 
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N, 
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300 

(76, 

69) 

320 

(75, 

69) 

310 

(74, 

69) 

340 

(73, 

69) 

290 

(78,70) 
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Table  3 

Position  of  Two  Blast  Sources  and  Seven  Microphones 


Blast  1 

84.68  Km, 

77. 48Km, 

291.10  m 

Microphonel 

84.68 

77.18 

289.74 

Microphone2 

82.20 

77.46 

271.88 

Microphone3 

79.25 

76.60 

303.95 

Microphone4 

76.78 

75.46 

333.28 

Microphone5 

74.54 

74.67 

336.47 

Microphone6 

71.95 

74.58 

320.83 

Microphone7 

69.81 

73.71 

330.22 

Blast  2 

69.28 

73.83 

333.21 
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Table  4 


Calculated  SELs  Generated  by  Blast  1  at  Fort  Leonard  Wood 


Blast  1,  generated  by  at  84.68  km,  77.48  km 

Window:  north:  78500.0  m 

south.  72500.0  m 

east:  86500.0  m 

west:  68500.0  m 

grid  size:  1000.00  m 

rows:  6 

cols:  1 8 


No- weighting,  1-lb  charge 


700 

749849 

709 

702 

890 

790 

980 

844 

1040  1007 

1036 

957 

858 

1073 

1048 

1078 

1048 

696 

605807 

621 

744 

590 

847 

737 

924 

902 

931 

694 

807 

1070 

1078 

1078 

1627 

1078 

759 

620807 

758 

705 

832 

967 

975 

830 

992 

908 

919 

778 

1000 

1032 

1048 

1078 

1048 

858 

703865 

838 

717 

426 

716 

988 

800 

1011 

1018 

1025 

974 

1017 

823 

856 

1079 

1067 

587 

802822 

871 

675 

778 

811 

991 

840 

1021 

889 

1017 

997 

975 

638 

878 

889 

853 

701 

792924 

924 

826 

868 

789 

888 

831 

954 

879 

976 

979 

727 

864 

680 

949 

915 

No-weighting,  5-lb  charge 


792 

829920 

782 

781 

974 

873  1063 

931 

1124  1084 

1119 

1036 

949 

1157 

1130 

1160 

1130 

784 

699888 

711 

840 

694 

915  810 

1000 

984  1016 

768 

890 

1153 

1152 

1160 

1705 

1160 

848 

718890 

833 

798 

908 

1051  1059 

906 

1069  1009 

1006 

863 

1067 

1109 

1130 

1160 

1130 

942 

782946 

944 

818 

494 

778  1069 

884 

1082  1084 

1096 

1054 

1092 

897 

939 

1153 

1142 

665 

893923 

944 

755 

853 

885  1073 

926 

1107  984 

1098 

1074 

1052 

712 

957 

974 

945 

806 

88TO11 

1019 

906 

963 

889  985 

918 

1025  961 

1051 

1053 

831 

924 

754 

1032 

1005 

c-weighting,  1-lb  charge 


673 

728830 

689 

679 

867 

767 

957 

819 

1017 

988 

1015 

936 

831 

1050 

1027 

1057 

1027 

671 

576786 

594 

714 

556 

830 

719 

906 

881 

908 

676 

785 

1048 

1060 

1057 

1610 

1057 

735 

589784 

743 

682 

809 

944 

953 

811 

973 

874 

895 

753 

985 

1012 

1027 

1057 

1027 

835 

682844 

801 

683 

410 

702 

968 

779 

993  1003 

1007 

954 

999 

805 

833 

1061 

1049 

566 

775789 

852 

652 

759 

789 

968 

817 

997 

861 

996 

977 

956 

620 

858 

865 

826 

667 

766902 

896 

806 

838 

758 

858 

808 

937 

856 

957 

961 

693 

852 

662 

926 

889 

c-weighting,  5-lb  charge 


752 

790888 

752 

744 

936 

828 

1024 

890 

1085 

1051 

1081 

1001 

903 

1119 

1092 

1122 

1092 

741 

651852 

666 

789 

639 

888 

781 

967 

946 

976 

737 

851 

1114 

1121 

1122 

1673 

1122 

805 

667852 

803 

757 

875 

1011 

1020 

873 

1036 

952 

964 

820 

1043 

1075 

1092 

1122 

1092 

903 

747910 

886 

765 

468 

756 

1032 

845 

1053 

1061 

1067 

1018 

1060 

866 

900 

1122 

1110 

630 

851878 

911 

716 

820 

850 

1035 

884 

1066 

934 

1060 

1040 

1018 

681 

922 

933 

899 

757 

844970 

973 

860 

913 

833 

935 

875 

995 

923 

1019 

1022 

774 

906 

723 

993 

960 
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Table  5 


Calculated  SELs  Generated  by  Blast  2  at  Fort  Leonard  Wood 


window:  north:  78500.0  m 

south:  72500.0  m 

east:  86500.0  m 

west:  68500.0  m 

grid  size:  1000.0  m 

rows:  6 

cols:  18 


No-weighting,  1 -pound  charge 


871 

952 

500  1039 

907 

700 

909 

701 

901 

652 

709 

768 

677 

680 

704 

718 

698 

704 

888 

889 

878 

817 

922 

717 

922 

567 

826 

500 

500 

726 

500 

765 

500 

500 

711 

500 

963 

1063 

1029 

850 

927 

883 

864 

574 

728 

703 

500 

749 

657 

657 

656 

720 

724 

681 

1018 

1048 

1027 

832 

824 

698 

816 

765 

545 

624 

500 

524 

500 

657 

500 

496 

707 

500 

1627 

1078 

830 

794 

833 

808 

708 

836 

671 

806 

622 

741 

500 

500 

500 

710 

775 

808 

1018 

1048 

817 

810 

823 

704 

783 

770 

692 

742 

734 

651 

738 

656 

632 

624 

726 

661 

No-weighting,  5-pound  charge 


948 

1033 

500  1123 

996 

786 

999 

794 

989 

749 

803 

858 

777 

775 

799 

79>i 

791 

799 

967 

972 

961 

911 

1011 

808 

1015 

682 

916 

500 

500 

828 

500 

865 

500 

500 

814 

500 

1033 

1138 

1104 

941 

1016 

967 

949 

676 

821 

797 

500 

846 

745 

744 

741 

814 

808 

784 

1100 

1130 

1106 

911 

916 

793 

911 

845 

654 

734 

500 

610 

500 

751 

500 

583 

787 

500 

1705 

1160 

901 

887 

920 

896 

793 

926 

763 

896 

721 

832 

500 

500 

500 

801 

863 

897 

1100 

1130 

887 

899 

909 

793 

875 

851 

777 

829 

813 

738 

838 

752 

736 

707 

834 

757 

C-weighting,  1 -pound  charge 


853 

930 

500 

1017 

882 

675 

883 

677 

876 

622 

681 

742 

644 

651 

674 

699 

671 

674 

686 

866 

855 

789 

897 

688 

893 

525 

800 

500 

500 

693 

500 

732 

500 

500 

679 

500 

948 

1045 

1010 

824 

901 

860 

840 

549 

700 

675 

500 

720 

632 

632 

633 

694 

703 

650 

997 

1027 

1007 

812 

797 

670 

788 

743 

509 

587 

500 

500 

500 

629 

500 

471 

686 

500 

1610 

1057 

813 

766 

809 

783 

684 

809 

643 

781 

590 

715 

500 

500 

500 

683 

750 

782 

997 

1027 

800 

785 

799 

678 

758 

749 

666 

718 

715 

626 

707 

625 

598 

603 

685 

632 

C-weighting,  5 

-pound 

charge 

915 

996 

500 

1983 

952 

740 

955 

745 

946 

700 

756 

810 

723 

728 

751 

760 

748 

748 

931 

933 

922 

863 

968 

764 

966 

617 

871 

500 

500 

770 

500 

817 

500 

500 

763 

500 

1006 

1106 

1072 

896 

972 

928 

908 

623 

774 

749 

500 

796 

702 

702 

700 

765 

768 

732 

1062 

1092 

1071 

876 

870 

744 

862 

809 

594 

672 

500 

569 

500 

703 

500 

541 

749 

500 

1673 

1122 

873 

840 

878 

853 

753 

881 

717 

852 

669 

787 

500 

500 

500 

755 

820 

853 

1062 

1092 

860 

856 

868 

749 

834 

812 

735 

786 

779 

696 

788 

705 

682 

666 

779 

712 

45 


47 
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7  SUMMARY 


This  report  describes  a  GTD  propagation  model  that  approximates  acoustic  wave  propagation  over 
hills  and  realistic  terrain.  The  realistic  terrain  data  can  be  either  input  by  the  user  or  supplied  by  a 
GRASS  geographic  information  system  data  base.  The  calculation  is  basically  a  point-to-point  calculation. 
When  the  terrain  data  are  input  by  the  user  and  one  observation  point  is  considered,  the  computation  is 
straightforward  and  very  quick.  The  user  can  obtain  not  only  the  total  sound  level,  but  also  the 
components  of  the  sound  levels  generated  by  different  mechanisms,  e.g.,  directly  propagated  from  source 
to  the  receiver,  reflected  from  the  ground  surface;  diffracted  from  hill  edges  and  so  on. 

The  GTD  model  was  interfaced  with  GRASS  programs  to  predict  the  SEL  for  an  explosive  source 
above  a  specified  terrain.  The  modified  system  calculates  SEL's  for  the  defined  area  to  produce  an  SEL 
map.  Sample  runs  on  data  from  the  Pinyon  Canyon,  CO  and  the  Fort  Leonard  Wood.  MO  regions  were 
described.  The  final  version  of  the  program  was  designed  for  simplicity.  The  user  supplies  a  source 
position,  charge  size,  and  an  expected  weighting,  and  then  the  computer  calculates  the  SEL  for  each  grid 
cell  of  a  given  terrain  area,  and  creates  a  contour  map  for  the  SEL  in  this  area.  For  this  study,  GRASS 
programs  were  altered  so  that  the  user  need  only  type  “gtdcom"  to  begin  the  program,  and  then  follow 
a  series  of  screen  prompts  to  enter  a  noise-weighting  index  (A  to  E)  and  an  effective  charge  weight  in 
pounds.  As  the  program  calculates,  the  SEL  results  for  each  cell  of  the  terrain  are  displayed  on  the  screen, 
so  that  the  user  may  monitor  the  process.  After  all  cell  calculations  are  completed,  a  hard  copy  of  the 
contour  map  of  SEL  values  is  prepared  and  sent  to  a  printer. 

Comparison  of  the  calculated  SEL  values  and  the  contour  maps  with  the  observed  data  shows  a 
tendency  for  agreement  between  the  variation  of  observed  SEL  with  distance,  and  (theoretical)  SEL 
predictions. 

The  program  calculations  used  in  this  study  were  based  on  ray  acoustic  theory.  The  application  of 
ray  theory  to  sound  waves  is  an  approximation,  especially  at  very  low  frequencies.  For  this  reason,  high 
accuracy  is  of  little  importance  for  the  other  aspects  of  the  calculation.  For  example,  the  diffraction  theory 
for  a  infinite  edge  was  used  in  the  situation  with  a  finite  edge  for  three-dimensional  terrain.  At  most, 
doubly-diffracted  rays  were  included;  triple  diffraction  rays,  and  so  on  were  excluded.  Although  we  found 
it  feasible  to  apply  geometric  ray  theory  to  sound  wave  propagation  over  realistic  terrain,  much 
improvement  would  be  expected  if  the  calculation  were  based  on  wave  theory  instead  of  the  approxima¬ 
tions  inherent  in  geometric  ray  theory. 

The  original  gtd  program  required  that  terrain  data  be  supplied  in  such  a  way  that  consecutive  pairs 
of  points  formed  parallel  lines.  This  requirement  set  a  strong  constraint  on  the  calculation,  with 
unavoidable  implications.  The  present  procedure  calculates  diffraction  and  reflection  simultaneously.  This 
procedure  may  not  be  best  when  applying  geometric  ray  theory  in  three  dimensions.  A  possible  alternative 
might  be  first  to  find  the  hills  (edges)  and  planes  for  the  terrain  area,  then  to  calculate  diffraction  rays  and 
reflected  rays,  and  finally  to  sum  them  up  to  get  SEL  values.  The  proposed  alternate  method  may  save 
computer  time  and  avoid  losing  rays  frxvn  any  plates  and  edges  in  the  considered  area.  One  way  to 
achieve  this  would  be  to  carefully  select  diffraction  edges  from  a  given  terrain  to  reduce  the  number  of 
edges  considered  and  to  include  all  possible  rays.  Only  important  edges  would  be  used  as  diffraction 
edges,  rather  than  every  cell  boundary  under  the  SOL  in  the  map. 
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This  study  revealed  several  sources  of  error  in  calculation  that  are  expected  to  be  reduced  in  future 

work: 


1.  To  save  computer  time  and  maintain  accuracy,  only  the  terrain  effects  under  the  direct  path  from 
the  source  to  the  observation  point  were  considered.  This  study  ignored  the  terrain  effects  in  the  backward 
direction  from  the  source-observation  line.  If  the  source  and  observation  lie  in  a  valley,  the  reflection  ray 
from  the  side  plates  may  be  significant.  It  would  not  be  very  difficult  to  include  the  effects  from  the  side 
plates  and  the  back  plates,  but  this  was  not  done  in  this  development. 

2.  The  sound  waves  produced  by  a  heavy  weapon  may  propagate  over  several  kilometers.  To  find 
out  the  SEL  for  the  whole  affected  region,  the  resolution  of  terrain  cells  should  not  be  so  fine  that  the 
number  of  cells  is  consequently  very  large.  To  model  the  effects  of  fine  structure  near  the  source  on  the 
noise  level,  the  cell  size  should  be  set  relatively  small  to  obtain  the  SEL  map  for  a  small  area.  However, 
the  accuracy  of  SEL  results  for  a  finer  resolution  is  not  necessarily  higher  than  a  rough  resolution. 
Furthermore,  the  resolution  should  be  set  in  relation  to  the  region  of  terrain  area.  The  gtd  programs 
demand  the  same  resolution  for  both  x  and  y  coordinates,  while  GRASS  may  supply  different  resolutions 
for  x  and  y  coordinates,  i.e.,  the  east-west  and  north-south  directions.  To  avoid  this  possible  discrepancy, 
the  user  should  ensure  that  the  resolution  and  the  sizes  of  the  terrain  area  in  both  x  and  y  directions  are 
the  same. 
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APPENDIX:  Program  Listing 


Subroutine  Page 

main.f*  main.w . 31 

atgn2.f  . 42 

blast.f** . 43 

blastl.f** . 44 

cot.f . 45 

csqr.f  . 46 

csx.f . 47 

cvn3d.f . 48 

dlogg.f  . 49 

dreal.f**  . 50 

f.f . 51 

fcrcfl.f  . 52 

image.f . 54 

P2intf**  . 56 

nwsrce.f  xx 

raydis.f . 58 

reOd.f . 59 

refvec.f . 60 

sddif.f*  sddif.w  . 61 

sgn.f . 65 

Si  .I**  . 66 

triple.! . 67 

Modified  subroutine. 


Newly  created  subroutine. 
Deleted  subroutine. 


Subroutine  Page 

tunder.f**  tunder.w  . 68 

v2dif/ . 82 

vamph3.f*  vamph3.w . 84 

vampha.f . 85 

vbeta.f* . 86 

vdrr.f  vdrr.w . 88 

vector.f . 91 

vengo.f  vengo.w  . 92 

vetdif.f  vetdif.w  . 99 

vexist.f .  110 

vgroua.f**  vgroua.w . 112 

vgroun.f* . 116 

vgtd.f . 120 

vlogil.f . 128 

vlogi2/*' . 129 

vmergh.f'*  . 130 

vrdd.f*  vrdd.w  . 131 

vrdr.f*  vrdr.w . 135 

vrrd.f  vrrd.w . 138 

vnT.f  vrrr.w . 141 

weight.f**  . 146 
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